Mendelian Randomization analysis of DeepSEA generated data

library(MendelianRandomization)

data_dir = "../../dat"
seq_predictions_all = read.csv(file.path(data_dir, "means_and_uncertainties_hopefully_correct.csv"))
seq_predictions = subset(seq_predictions_all, seq_num == 25)

str(seq_predictions)

'data.frame':   3001 obs. of  6 variables:
 $ seq_num    : int  25 25 25 25 25 25 25 25 25 25 ...
 $ mut        : Factor w/ 3000 levels "0A","0C","0G",..: 12 1 2 3 342 343 344 674 675 676 ...
 $ X_pred_mean: num  5.66e-06 2.34e-05 4.55e-05 -3.47e-05 1.45e-04 ...
 $ X_pred_var : num  3.17e-05 3.95e-05 6.99e-05 3.93e-05 9.92e-05 ...
 $ Y_pred_mean: num  4.68e-06 1.52e-05 8.47e-05 -9.28e-06 8.67e-05 ...
 $ Y_pred_var : num  2.52e-05 3.08e-05 8.28e-05 2.55e-05 9.15e-05 ...

bx <- unlist(seq_predictions["X_pred_mean"])
bxse <- unlist(seq_predictions["X_pred_var"])
by <- unlist(seq_predictions["Y_pred_mean"])
byse <- unlist(seq_predictions["Y_pred_var"])

MRInputObjectEx <- mr_input(bx = bx, 
                          bxse = bxse, 
                          by = by, 
                          byse = byse)

cis = c()
for (seq in (1:50)) {
  seq_predictions = subset(seq_predictions_all, seq_num == seq)
  
  bx <- unlist(seq_predictions["X_pred_mean"])
  bxse <- unlist(seq_predictions["X_pred_var"])
  by <- unlist(seq_predictions["Y_pred_mean"])
  byse <- unlist(seq_predictions["Y_pred_var"])
  MRInputObject <- mr_input(bx = bx, 
                          bxse = bxse, 
                          by = by, 
                          byse = byse)
  EggerObject <- mr_egger(MRInputObject, robust = FALSE, alpha = .05)
  # cis = append(cis, c(EggerObject$CILower.Est, EggerObject$CIUpper.Est))
  print(seq)
  print(c(EggerObject$CILower.Est, EggerObject$CIUpper.Est))
  print(c(EggerObject$I.sq))
}

[1] 1
[1] 0.05687267 0.06067051
[1] 0
[1] 2
[1] 0.3306068 0.3624528
[1] 0.2105841
[1] 3
[1] 0.03796090 0.04176579
[1] 0.3214789
[1] 4
[1] 0.1549294 0.1735800
[1] 0.2925369
[1] 5
[1] 0.1533373 0.1661913
[1] 0.07879359
[1] 6
[1] 0.1685509 0.1815121
[1] 0
[1] 7
[1] 0.2675980 0.2861198
[1] 0
[1] 8
[1] 0.2909661 0.3147951
[1] 0.1831951
[1] 9
[1] 0.1449050 0.1556574
[1] 0.1622787
[1] 10
[1] 0.1257318 0.1356557
[1] 0
[1] 11
[1] 0.1473885 0.1603717
[1] 0
[1] 12
[1] 0.04852230 0.05193937
[1] 0
[1] 13
[1] 0.2326184 0.2606935
[1] 0
[1] 14
[1] 0.4096944 0.4500918
[1] 0
[1] 15
[1] 0.1006563 0.1091880
[1] 0.02123769
[1] 16
[1] 0.03630137 0.03925981
[1] 0.2420155
[1] 17
[1] 0.08023490 0.08748782
[1] 0.1503298
[1] 18
[1] 0.2139596 0.2311665
[1] 0
[1] 19
[1] 0.2752876 0.2981095
[1] 0
[1] 20
[1] 0.2259239 0.2427902
[1] 0
[1] 21
[1] 0.01434114 0.01553469
[1] 0
[1] 22
[1] 0.1928645 0.2083450
[1] 0
[1] 23
[1] 0.05310580 0.05783676
[1] 0
[1] 24
[1] 0.3177343 0.3388181
[1] 0
[1] 25
[1] 0.2673437 0.2889935
[1] 0.2421144
[1] 26
[1] 0.02250383 0.02425483
[1] 0
[1] 27
[1] 0.05180822 0.05611104
[1] 0
[1] 28
[1] 0.05388898 0.05843682
[1] 0
[1] 29
[1] 0.04815811 0.05224462
[1] 0
[1] 30
[1] 0.08851500 0.09372713
[1] 0.174372
[1] 31
[1] 0.3338592 0.3653642
[1] 0
[1] 32
[1] 0.03867074 0.04214653
[1] 0
[1] 33
[1] 0.4801847 0.5207160
[1] 0
[1] 34
[1] 0.03362065 0.03616059
[1] 0
[1] 35
[1] 0.05669501 0.06187179
[1] 0
[1] 36
[1] 0.1737785 0.1876079
[1] 0.1640396
[1] 37
[1] 0.08575466 0.09214800
[1] 0
[1] 38
[1] 0.06662938 0.07251736
[1] 0
[1] 39
[1] 0.1864075 0.2002264
[1] 0
[1] 40
[1] 0.2280670 0.2459672
[1] 0
[1] 41
[1] 0.02417243 0.02568200
[1] 0
[1] 42
[1] 0.08946369 0.09645143
[1] 0
[1] 43
[1] 0.1533693 0.1655198
[1] 0
[1] 44
[1] 0.04358466 0.04715561
[1] 0
[1] 45
[1] 0.1027557 0.1105777
[1] 0
[1] 46
[1] 0.1188304 0.1293053
[1] 0
[1] 47
[1] 0.02712174 0.02866630
[1] 0.05283552
[1] 48
[1] 0.1388547 0.1522203
[1] 0
[1] 49
[1] 0.4762848 0.5147806
[1] 0
[1] 50
[1] 0.09995734 0.10829652
[1] 0.09234253

# dim(cis) <- c(2, 50)
# print(cis)

mr_plot(MRInputObjectEx)

install.packages("devtools")

Installing package into ‘/home/stephenmalina/R/x86_64-pc-linux-gnu-library/3.6’
(as ‘lib’ is unspecified)
also installing the dependency ‘remotes’

trying URL 'https://cloud.r-project.org/src/contrib/remotes_2.1.1.tar.gz'
Content type 'application/x-gzip' length 134042 bytes (130 KB)
==================================================
downloaded 130 KB

trying URL 'https://cloud.r-project.org/src/contrib/devtools_2.2.2.tar.gz'
Content type 'application/x-gzip' length 375464 bytes (366 KB)
==================================================
downloaded 366 KB

* installing *source* package ‘remotes’ ...
** package ‘remotes’ successfully unpacked and MD5 sums checked
** using staged installation
** R
** inst
** byte-compile and prepare package for lazy loading
** help
*** installing help indices
** building package indices
** installing vignettes
** testing if installed package can be loaded from temporary location
** testing if installed package can be loaded from final location
** testing if installed package keeps a record of temporary installation path
* DONE (remotes)
* installing *source* package ‘devtools’ ...
** package ‘devtools’ successfully unpacked and MD5 sums checked
** using staged installation
** R
** inst
** byte-compile and prepare package for lazy loading
** help
*** installing help indices
*** copying figures
** building package indices
** installing vignettes
** testing if installed package can be loaded from temporary location
** testing if installed package can be loaded from final location
** testing if installed package keeps a record of temporary installation path
* DONE (devtools)

The downloaded source packages are in
    ‘/tmp/RtmpfKGeTB/downloaded_packages’

library(devtools)

Loading required package: usethis

install_github("WSpiller/RadialMR")

Skipping install of 'RadialMR' from a github remote, the SHA1 (18a7a6e9) has not changed since last install.
  Use `force = TRUE` to force installation

plot(by/bx, bx/byse, xlim=c(min((by-2*byse)/bx), max((bx+2*byse)/bx)), ylim=c(0, max(bx/byse)))
for (j in 1:length(bx)) {
lines(c((by[j]-1.96*byse[j])/bx[j], (by[j]+1.96*byse[j])/bx[j]),
c(bx[j]/byse[j], bx[j]/byse[j]))
}

abline(h=0, lwd=1); abline(v=0, lwd=1)

library(RadialMR)
r_input <- format_radial(bx, by, bxse, byse)
funnel_radial(egger_radial(r_input, alpha = .05))

funnel_radial(ivw_radial(r_input, alpha = .05))

Weights not specified: Adopting modified second-order weights

Radial IVW

Residual standard error: 0.948 on 3000 degrees of freedom

F-statistic: 2820.08 on 1 and 3000 DF, p-value: <2e-16
Q-Statistic for heterogeneity: 2698.411 on 3000 DF , p-value: 0.9999713

 Outliers detected 
Number of iterations = 6

Error: Tibble columns must have consistent lengths, only values of length one are recycled:
* Length 4: Columns `x`, `xend`
* Length 3001: Column `PANEL`

plot_radial(egger_radial(r_input, alpha = .05))